
## ---------------------- sim1 and sim2 block diagnoal w matrix

rm(list = ls())
## compile functions 
library(bpNet)
##library(mcmcr)
library(coda)
library(MASS)
library(abind)

## Set working directory where the data, code, and output are stored
setwd('~/Desktop/PSRMreplicationLP/')

set.seed(123456789)

## ------------------------- ##
source("code/simulated20529.R")   ## code for generating simulated data
source("code/net_plot.R")

## the correlation between factor loadings and the network ties
## we generate three datasets with strength of confounding 0.3, 0.6, and 0.9, respectively

N <- 200 ## number of units
Ng <- 20 ## number of groups 
TT <- 50 ## time period
unit <- N / Ng

## 1. gen W matrix 
w0 <- matrix(0, unit, unit)
diag(w0) <- 1


w <- (matrix(1, unit, unit) - w0)/(unit - 1)
w1 <- matrix(0, N, N)
for (i in 1:10) {
	w1[((i-1)*unit + 1):(i*unit),((i-1)*unit + 1):(i*unit)] <- w
}


W <- array(0, dim = c(N, N, TT))
for(i in 1:TT) {
    W[,,i] <- w1
}

save(W, file="data/SimData/blockW.RData")
## 2. gen time-varying covar that explains rho_t
rhoZ <- cbind(1, rnorm(TT))


## 3. gen other parameters

## covar
p <- 2
beta <- as.matrix(c(5, 1, 3, -2, 2))

## rho_t
rhoA <- as.matrix(c(0.3, -0.2))
kappa <- 0.3
sigma_n2 <- 0.36

## factor numbers
r <- 2

## error variance
sigma2 <- 1

## contextual effects
contextual.effect <- "both"
force <- "both"

## Tie-Homophily 0.3

    ## generate a dateset for Study I: varying rho_t 
    sim <- spatial.simulate(W = W, ## W matrix
                        beta = beta,
		                p = p, 
		                r = r,
		                contextual.effect = contextual.effect,
		                force = force,
		                sigma2 = sigma2,
		                kappa = kappa, 
		                rhoZ = rhoZ,
		                rhoA = rhoA, 
		                sigma_n2 = sigma_n2, corH=0.3)
    ## Data for simulation study I
    sim <- list(sim, rhoZ)
    save(sim, file="data/SimData/sim1cor3.RData")

   ## Dataset for Study II: rho_t=0
   sim2 <- spatial.simulate(W = W, ## W matrix
                        beta = beta,
		                p = p, 
		                r = r,
		                contextual.effect = contextual.effect,
		                force = force,
		                sigma2 = sigma2,
		                kappa = 0, 
		                rhoZ = rhoZ,
		                rhoA = as.matrix(c(0, 0)),  # no alpha_1 (the true value of alpha1 is zero)
		                sigma_n2 = 0, corH=0.3)
      ## Data for simulation study II
       sim2 <- list(sim2, rhoZ)
       save(sim2, file="data/SimData/sim2cor3.RData")

rm(list = ls())
## compile functions 
library(bpNet)
##library(mcmcr)
library(coda)
library(MASS)
library(abind)

## Set working directory where the data, code, and output are stored
## setwd('~/Desktop/PSRMreplicationLP/code')

set.seed(123456789)

## ------------------------- ##
source("code/simulated20529.R")   ## code for generating simulated data
source("code/net_plot.R")

## the correlation between factor loadings and the network ties
## we generate three datasets with strength of confounding 0.3, 0.6, and 0.9, respectively

N <- 200 ## number of units
Ng <- 20 ## number of groups 
TT <- 50 ## time period
unit <- N / Ng

## 1. gen W matrix 
w0 <- matrix(0, unit, unit)
diag(w0) <- 1


w <- (matrix(1, unit, unit) - w0)/(unit - 1)
w1 <- matrix(0, N, N)
for (i in 1:10) {
	w1[((i-1)*unit + 1):(i*unit),((i-1)*unit + 1):(i*unit)] <- w
}


W <- array(0, dim = c(N, N, TT))
for(i in 1:TT) {
    W[,,i] <- w1
}

## 2. gen time-varying covar that explains rho_t
rhoZ <- cbind(1, rnorm(TT))


## 3. gen other parameters

## covar
p <- 2
beta <- as.matrix(c(5, 1, 3, -2, 2))

## rho_t
rhoA <- as.matrix(c(0.3, -0.2))
kappa <- 0.3
sigma_n2 <- 0.36

## factor numbers
r <- 2

## error variance
sigma2 <- 1

## contextual effects
contextual.effect <- "both"
force <- "both"
## Tie-Homophily 0.6

    ## generate a dateset for Study I: varying rho_t 
    sim <- spatial.simulate(W = W, ## W matrix
                        beta = beta,
		                p = p, 
		                r = r,
		                contextual.effect = contextual.effect,
		                force = force,
		                sigma2 = sigma2,
		                kappa = kappa, 
		                rhoZ = rhoZ,
		                rhoA = rhoA, 
		                sigma_n2 = sigma_n2, corH=0.6)
    ## Data for simulation study I
    sim <- list(sim, rhoZ)
    save(sim, file="data/SimData/sim1cor6.RData")

   ## Dataset for Study II: rho_t=0
   sim2 <- spatial.simulate(W = W, ## W matrix
                        beta = beta,
		                p = p, 
		                r = r,
		                contextual.effect = contextual.effect,
		                force = force,
		                sigma2 = sigma2,
		                kappa = 0, 
		                rhoZ = rhoZ,
		                rhoA = as.matrix(c(0, 0)),  # no alpha_1 (the true value of alpha1 is zero)
		                sigma_n2 = 0, corH=0.6)
      ## Data for simulation study II
       sim2 <- list(sim2, rhoZ)
       save(sim2, file="data/SimData/sim2cor6.RData")

## Clean the environment to avoid any correlation between two datasets

rm(list = ls())
## compile functions 
library(bpNet)
##library(mcmcr)
library(coda)
library(MASS)
library(abind)

## Set working directory where the data, code, and output are stored
set.seed(123456789)

## ------------------------- ##
source("code/simulated20529.R")   ## code for generating simulated data
source("code/net_plot.R")

## the correlation between factor loadings and the network ties
## we generate three datasets with strength of confounding 0.3, 0.6, and 0.9, respectively

N <- 200 ## number of units
Ng <- 20 ## number of groups 
TT <- 50 ## time period
unit <- N / Ng

## 1. gen W matrix 
w0 <- matrix(0, unit, unit)
diag(w0) <- 1


w <- (matrix(1, unit, unit) - w0)/(unit - 1)
w1 <- matrix(0, N, N)
for (i in 1:10) {
	w1[((i-1)*unit + 1):(i*unit),((i-1)*unit + 1):(i*unit)] <- w
}


W <- array(0, dim = c(N, N, TT))
for(i in 1:TT) {
    W[,,i] <- w1
}

## 2. gen time-varying covar that explains rho_t
rhoZ <- cbind(1, rnorm(TT))


## 3. gen other parameters

## covar
p <- 2
beta <- as.matrix(c(5, 1, 3, -2, 2))

## rho_t
rhoA <- as.matrix(c(0.3, -0.2))
kappa <- 0.3
sigma_n2 <- 0.36

## factor numbers
r <- 2

## error variance
sigma2 <- 1

## contextual effects
contextual.effect <- "both"
force <- "both"

## Clean the environment to avoid any correlation between two datasets

## Tie-Homophily 0.9

    ## generate a dateset for Study I: varying rho_t 
    sim <- spatial.simulate(W = W, ## W matrix
                        beta = beta,
		                p = p, 
		                r = r,
		                contextual.effect = contextual.effect,
		                force = force,
		                sigma2 = sigma2,
		                kappa = kappa, 
		                rhoZ = rhoZ,
		                rhoA = rhoA, 
		                sigma_n2 = sigma_n2, corH=0.9)
    ## Data for simulation study I
    sim <- list(sim, rhoZ)
    save(sim, file="data/SimData/sim1cor9.RData")

   ## Dataset for Study II: rho_t=0
   sim2 <- spatial.simulate(W = W, ## W matrix
                        beta = beta,
		                p = p, 
		                r = r,
		                contextual.effect = contextual.effect,
		                force = force,
		                sigma2 = sigma2,
		                kappa = 0, 
		                rhoZ = rhoZ,
		                rhoA = as.matrix(c(0, 0)),  # no alpha_1 (the true value of alpha1 is zero)
		                sigma_n2 = 0, corH=0.9)
      ## Data for simulation study II
       sim2 <- list(sim2, rhoZ)
       save(sim2, file="data/SimData/sim2cor9.RData")


##------- Study III:  Data generation  ---------##
rm(list = ls())
## compile functions 
library(bpNet)
#library(mcmcr)
library(mvtnorm)
library(coda)
set.seed(123456789)

## ------------------------- ##
source("code/simulated30529.R")
source("code/net_plot.R")

## Load the network based on cleaned ICEWS 
load("data/W.RData")
## ------------- gen W matrix ------------ ## 

## 1. get the first 50 periods
W <- Y_ICEWS[,,1:50]
N <- dim(W)[1]
TT <- dim(W)[3]

rmn <- 3 ## remove 3 units that connect with almost every unit 
W2 <- array(NA, dim = c(N-rmn, N-rmn, TT))

for (i in 1:TT) {
	subw <- W[,,i]
	subw <- subw[-c(47, 8, 34), ]
	subw <- subw[, -c(47, 8, 34)]

	## only keep 3 units that have largest weights as neighbors
        ## and standardize the rows
	for (j in 1:(N-rmn)) {

		subv <- cbind(1:(N-rmn), subw[j, ])
		subv <- subv[order(subv[, 2]), ]

		pos1 <- subv[1:(N-rmn-3), 1]
		## pos2 <- subv[(N-4):N, 1]

		subw[j, pos1] <- 0
		## subw[j, pos2] <- 1

		subw[j,] <- subw[j,] / sum(subw[j,])
	}
	W2[,,i] <- subw
}

## obstain the w matrix for use 
W <- W2

save(W, file="data/SimData/icewsW.RData")

## generate time-varying covar that explains rho_t
rhoZ <- cbind(1, rnorm(TT))

## ------------- set other parameters ------------ ## 

## covar
p <- 2
beta <- as.matrix(c(5, 1, 3, -2, 2))

## rho_t
rhoA <- as.matrix(c(0.3, -0.2))
kappa <- 0.3
sigma_n2 <- 0.36

## factor numbers
r <- 2

## error variance
sigma2 <- 1

## contextual effects
contextual.effect <- "both"
force <- "both"

## Generate simulation data
sim3 <- spatial.simulate(W = W, 
                        beta = beta,
		                p = p, 
		                r = r,
		                contextual.effect = contextual.effect,
		                force = force,
		                sigma2 = sigma2,
		                kappa = kappa, 
		                rhoZ = rhoZ,
		                rhoA = rhoA, 
		                sigma_n2 = sigma_n2)

  sim3 <- list(sim3, rhoZ)
  save(sim3, file="data/SimData/sim3.RData")


